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Abstract 

Improvements to a time-accurate approximate 
factorization (AF) algorithm have been implemented for 
steady and unsteady transonic analysis of realistic aircraft 
configurations. These algorithm improvements have been 
made to the CAP-TSD (Computational Aeroelasticity 
Erogram * Xransonic Email Disturbance) code developed 
recently at NASA Langley Research Center. The code permits 
the aeroetastic analysis of complete aircraft in the flutter 
critical transonic speed range. The AF algorithm of the CAP* 
TSD code solves the unsteady transonic small-disturbance 
equation. The algorithm improvements include: an 

Engquist-Osher (E-O) type-dependent switch to more 
accurately and efficiently treat regions of supersonic flow, 
extension of the E-0 switch for second-order spatial 
accuracy in these regions, nonreflecting far field boundary 
conditions for more accurate unsteady applications, and 
several modifications which accelerate convergence to 
steady-state. Calculations are presented for several 
configurations including the General Dynamics one-ninth 
scale F-16C aircraft model to evaluate the algorithm 
modifications. The modifications have significantly 
improved the stability of the AF algorithm and hence the 
reliability of the CAP-TSD code in general. The paper 
presents detailed descriptions of the algorithm 
improvements along with results and comparisons which 
demonstrate the improved stability, accuracy, and efficiency 
of the CAP-TSD oode. 

Nomencfatme 

c airfoil chord 

Cia unsteady lift-curve slope 

Cr wing reference chord 

Cp pressure coefficient 

k reduced frequency. <ocr/2U 

M freestream Mach number 

NSUP number of supersonic points 

R residual 

t time, nondimensionalized by freestream speed and 

wing reference chord 
U freestream speed 

a instantaneous angle of attack 

oto mean angle of attack 

at amplitude of pitch oscillation 

y ratio of specific heats 

At nondimensional time step 

ii fractional semispan 

+ disturbance velocity potential 

Co angular frequency 

Subscript* 

t tail 

w wing 


Introduction 

Presently, considerable research is being conducted to 
develop finite-difference computer codes for calculating 
transonic unsteady aerodynamics for aeroetastic 
applications. 1 These computer codes are being developed to 
provide accurate methods of calculating unsteady airloads for 
the prediction of aeroetastic phenomena such as flutter and 
divergence. For example, the CAP-TSD 2 unsteady transonic 
small-disturbance (TSD) code was recently developed for 
transonic aeroelastic analyses of complete aircraft 
configurations. The name CAP-TSD is an acronym for 
Computational Aeroelasticity Erogram - Xransonic Email 
Disturbance. The new code permits the calculation of 
unsteady flows about complete aircraft for aeroelastic 
analysis in the flutter critical transonic speed range. The 
code can treat configurations with arbitrary combinations of 
lifting surfaces and bodies including canard, wing, tail, 
control surfaces, tip launchers, pylons, fuselage, stores, and 
nacelles. In Ref. 2, steady and unsteady pressures were 
presented for several complex aircraft configurations which 
demonstrated the geometrical applicability of CAP-TSD. 
These calculated results were in good agreement with 
available experimental pressure data which validated CAP- 
TSD for multiple component applications with mutual 
aerodynamic interference effects. Preliminary aeroelastic 
applications of CAP-TSD were presented in Ref. 3 for a 
simple well-defined wing case. The case was selected as a 
first step toward performing aeroelastic analyses for 
complete aircraft configurations. The calculated flutter 
boundaries compared well with the experimental data for 
subsonic as well as supersonic freestream Mach numbers, 
which gives confidence in CAP-TSD for aeroelastic 
prediction. 

The CAP-TSD code uses a time accurate approximate 
factorization (AF) algorithm recently developed by Batina 4 
for solution of the unsteady TSD equation. The AF algorithm 
involves a Newton linearization procedure coupled with an 
internal iteration technique. In Ref. 4, the algorithm was 
shown to be efficient for application to steady or unsteady 
transonic flow problems. It can provide accurate solutions 
in only several hundred time steps, yielding a significant 
computational cost savings when compared to alternative 
metnods. For reasons of practicality and affordability, an 
efficient algorithm and a fast computer code are 
requirements for realistic aircraft applications. 

The purpose of this paper is to describe recent changes to 
the CAP-TSD code which have significantly improved the 
stability of the AF algorithm and the accuracy of the results. 
The algorithm modifications include: (1) improved type- 

dependent differencing to treat regions of supersonic flow, 
(2) extension of the type -dependent differencing for second- 
order spatial accuracy, (3) nonreflecting far field boundary 
conditions for unsteady applications, and (4) several 
modifications to accelerate convergence to steady state. The 
paper presents detailed descriptions of these algorithm 
improvements along with results and comparisons which 
assess the improved stability, accuracy, and efficiency of the 
CAP-TSD oode. 


Algorithm .Improvements 


Stwdy-SUtfl- fiflmtoiflftncft Acceleration 


Enaquist-Os her Type-Dependent Switch 

Algorithms based on the TSO equation typically use 
central differencing in regions of subsonic flow and upwind 
differencing in regions of supersonic flow. This, of course, 
allows for the correct numerical description of the physical 
domain of dependence. The original CAP-TSD code of Ref. 2 
used the Murman5 type-dependent switch to change the 
spatial differencing. The Murman switch, however, admits 
nonphysical expansion shocks as a part of the solution and 
has been shown to be less stable than monotone methods. 6 * 7 
For example, unsteady results for a NACA 64A006 airfoil 
were presented in Ref. 7 which demonstrated an order of 
magnitude increase in time step using a monotone algorithm. 
Therefore, an Engquist-Osher (E-O) monotone switch, 
similar to that of Ref. 6, has been incorporated within the AF 
algorithm of the CAP-TSD code. The E-0 switch is based on 
sonic reference conditions and does not admit expansion 
shocks as part of the solution. Use of the E-0 switch also 
generally increases computational efficiency because of the 
larger time steps which may be taken. Mathematical details 
of the required algorithm changes are described in a 
subsequent section. 


Second-Order Accurate Soailal .-Differencing 

Most TSD algorithms are only first-order-accurate 
spatially in regions of supersonic flow. This is due to the 
first-order upwind differencing that is typically used to 
treat these regions. Use of second-order upwind differencing 
has been shown to improve the accuracy of the solution while 
retaining the numerical stability of the first order 
method. 0 Consequently, the E O type-dependent switch of 
the AF algorithm has been extended for second-order spatial 
accuracy in supersonic regions of the flow. Comparisons of 
results obtained using first-order and second order 
differencing, to be presented, demonstrate the improved 
accuracy of the second order method. 


Finally, several algorithm changes have been made to 
accelerate convergence to steady-state. Besides the E-0 
switch, these changes include: (1) deletion of the time- 
dependent terms from the residual of the AF algorithm, (2) 
deletion of atl of the time-derivatives of the TSO equation, 
and (3) over-relaxation of the residual. The effects of each 
of these modifications on the steady-state convergence are 
demonstrated in the results presented herein. 


IranM.Qlc_SmflilrP)8tjurbancfl, Equation 

The flow is assumed to be governed by the general 
frequency modified TSO potential equation which may be 
written in conservation law form as 
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where 


f 0 = - A V B< \ 

(2a) 

f, = E0„ + F$; + G<f 

(2b) 

f 2 

(2c) 

f 3 = 0. 

(2d) 

The coefficients A. B, and E are defined as 


A = M 2 , B = 2M 2 , E = 1 - M 2 

(3) 


Several choices are available for the coefficients F, G, and H 
depending upon the assumptions used in deriving the TSD 
equation 2 The coefficients are herein defined as 


Nonreflectina Far Field Boundary Conditions 

For unsteady applications, the far field boundary 
conditions can have a significant influence on the accuracy of 
Ihe solution. Steady state boundary conditions are 
inadequate for unsteady calculations, since disturbances 
reaching the boundaries are reflected back into the 
computational domain. These reflected disturbances can 
propagate into the near field and thus produce inaccurate 
results. One solution to this problem is to locate the grid 
boundaries far away to minimize the effect of the boundary 
conditions. This is generally not an acceptable remedy 
because ol the higher computational cost which results from 
an increased number of grid points required to discretize a 
larger computational domain. The more appropriate solution 
is the use of nonreflecting far field boundary conditions 
which absorb most of the waves that are incident on the 
boundaries and consequently allow the use of smaller 
computational grids. 9 Nonretlecting boundary conditions 
similar to those of Whitlow 9 have been incorporated within 
the CAP-TSD code. These boundary conditions are consistent 
with the AF solution procedure and are described in more 
detail below. Results obtained with and without the 
nonreflecting boundary conditions are presented which 
demonstrate their effectiveness. 
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Approximate- factorization .Algorithm 

An approximate factorization algorithm was developed 4 
to solve Ihe modified TSD equation (Eq. (1)). In this section, 
the AF algorithm is described. 


Qgngrgl J kaP UP liPD 

The AF algorithm consists of a Newton linearization 
procedure coupled with an internal iteration technique. For 
unsteady flow calculations, the solution procedure involves 
two steps. First, a time linearization step (described below) 
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is performed to determine an estimate of the potential field. 

Second, internal iterations are performed to provide time < 

accurate modeling of the flow field. Specifically, the TSD F 2 « — (1 +Hl*V) (8e) 

equation (Eq. (1)) is written in general form as 5* * 


R(0 Bt ‘),O 


(5) 



(8f) 


where represents the unknown potential field at time 
level (n+1). The solution to Eq. (5) is then given by the 
Newton linearization of Eq. 5 about ♦* 


At 2 , A 20* - 50 n + V • <t>" 2 
X2A ^ 


R(40 + <— ) »A0«O (6) 

d0 ** 

In Eq. (6), p' is the currently available value of ♦ n+t and 
Ap «pn + i . During convergence of the iteration 
procedure, Ap will approach zero so that the solution will be 
given by p n *i * p*. In general, only one or two iterations 
are required to achieve acceptable convergence. For steady 
flow calculations, iterations are not used since time accuracy 
is not necessary when marching to steady state. 
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Mathematica l Formulation 

The AF algorithm is formulated by first approximating 
the time derivative terms (pit andPxt terms) by second- 
order accurate finite -difference foimulae. The TSO equation 
is rewritten by substituting f f and neglecting 

squares of derivatives of Ap which is equivalent to applying 
Eq. (6) term by term The resulting equation is then 
rearranged and the left-hand side is approximately factored 
into a triple product of operators yielding 


where 
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In Eq. (7) a is a relaxation parameter which is normally set 
equal to 1.0. To accelerate convergence to steady-state, the 
residual R may be over relaxed using o > 1 . Equation (7) is 
solved using three sweeps through the grid by sequentially 
applying the operators 1 4 , U), and U; as 


£ - sweep: 

Ap = - a R 

(9a) 

q - sweep. 
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(9b) 

£ - sweep: 
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Further details of the algorithm development and solution 
procedure may be found in Ref 4. 
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Enqqutst-Oiher Tvoa-D apandent Switch 

An Engquist-Osher type -dependent mixed difference 
operator has been implemented in the AF algorithm to treat 
supersonic regions of the flow. The E-0 switch is based on 
sonic reference conditions and is applied to both sides of Eq. 
(7). For example, in the residual (Eq. (8g|) the terms that 
are upwind biased at supersonic points are defined by 
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where 


^ V F ^ 2) Cl/2 lt/2 (10) 


f i-l/2“ Eu i-l/2 + FU i.|/2 


( 11 «) 


Similar modifications to the left-hand side of Eq. (7) result 
in a pentadiagonal system of equations for subsonic flows 
with embedded supersonic regions and a tridiagonaf system of 
equations for purely subsonic flows. Furthermore, the 
treatment of the $xi term in the TSD equation is only first- 
order accurate in space because of the one-sided differencing 
used. Similar to Ref. 8. the +xt term is backward differenced 
to enhance diagonal dominance and consequently maintain 
numerical stability. 


f i.|/2 = Eu i../2 + Fu i-l/2 
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Boundary Condition* 


a. f : 


' f i-l/2 


4 i ,/2 o/axtw't,.,) 


«- A 
A. f 


A A 

f i l/2 ' f i-3/2 


4 M/ * 


+ <U i i/2 - u > 


U M/2 = U+E i.l/2< U M/2- U > 


(He) 


(lid) 


(lie) 


(lit) 


Flow-tanpency. - The flow tangency boundary 
conditions are imposed along the mean plane of the respective 
lifting surfaces and the wakes are assumed to be planar 
extensions from the trailing edges to the downstream 
boundary of the finite-difference grid. The numerical 
implementation of these conditions? allows for copianar as 
well as non-coplanar combinations of horizontal (canard, 
wing, horizontal tail, launchers) and vertical (pylons, 
vertical tail) surfaces. Bodies such as the fuselage, stores, 
and nacelles are treated using simplified boundary conditions 
on a prismatic surface rather than on the true surface? The 
method is consistent with the small-disturbance 
approximation and treats bodies with sufficient accuracy to 
obtain the correct global effect on the flow field without the 
use of special grids or complicated coordinate 
transformations. 


U i.|/2 = S »*i 




u= sonic value of 


e i 1/2 = 1 if U i l/2 >U 


0 if U i l/2 — U 
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In Eqs. (11) the j and k subscripts corresponding to the 
spanwise and vertical directions, respectively, have been 
omitted for clarity. Similar differencing is used on the left- 
hand side of Eq. (7) where the first two terms of Ft (Eq. 8d) 
are upwind biased at supersonic points. 


Sacond-Ordar- Accurate Spatial BUffMtafl 

The AF algorithm with the E-0 switch as defined by Eq. 
(to) is only first-order accurate in supersonic regions ot 
the flow. To achieve second-order accuracy at supersonic as 
well as subsonic points, Eq. (10) is extended as 




Far Field. - The conditions imposed upon the outer 
boundary of the computational region are similar to the 
nonreflecting boundary conditions reported by Whitlow 9 
The conditions employed here are given by 


Upstream: 


Above: 


Below: 


Right spanwise: 


Left spanwise: 


Symmetry plane: 


♦ = <) 

(13a) 


(13b) 

fw° 

(13c) 


(13d) 

y "° 

(I3e) 

yV* y =0 

( 1 3 f ) 

(for full-span modeling) 

0 =0 

(I3g) 


♦ A *V‘-.«V‘t/2 (12) 


(for half-span modeling) 
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where C - E + 2F*x and D - v 4A + . These boundary 

conditions are numerically imposed by redefining the L&, Li), 
and L; operators in Eq. (7) as welt as the right-hand side R, 

at the appropriate grid points. The equation to be solved at 
boundary grid points may then be written symbolically as 

L t L L A<(> = - o R (14) 

where the "tilde" indicates that the quantity has been 
rewritten to account for the boundary conditions. For 


example, along the downstream boundary the three operators 
and right-hand side are defined as 

7 , 4 At -B D d 

L , = 1 +^-(7T+-=) \ t — 

4 3 c Vc 

(15a) 

L =1 
n 

(15b) 

ii 

(15c) 

R = i (30* - 40 - ^ 

(I5d) 


Time-Linearization Step 

An initial estimate of the potentials at time level (n+1) 
is required to start the iteration process. This estimate is 
provided by performing a time-linearization calculation. 
The equations governing the time-linearization step are 
derived in a similar fashion as the equations for iteration. 
The only difference is that the equations are formulated by 
linearizing about time level (n) rather than the iterate level 

n. 


CAP-TSP Code 

The AF algorithm has been used as the basis of the CAP- 
TSD code for transonic unsteady aerodynamic and aeroelastic 
analysis of realistic aircraft configurations. The code can 
treat configurations with arbitrary combinations of lifting 
surfaces and bodies including canard, wing, tail, control 
surfaces, tip launchers, pylons, fuselage, stores, and 
nacelles. The present capability has the option of half-span 
modeling (Eq (13g)) for symmetric cases or full-span 
modeling (Eq. ( 1 3f )) to allow the treatment of 
antisymmetric mode shapes, fuselage yaw, or unsymmetric 
configurations such as an oblique wing or unsymmetric wing 
stores. Steady and unsteady CAP-TSD pressures for several 
realistic aircraft configurations, including comparisons 
with experimental data, were presented in Ref. 2. The 
calculated results were in good agreement with the 
experimental pressure data which validated CAP-TSD for 
multiple component applications with mutual aerodynamic 
interference effects. Preliminary aeroelastic applications of 
CAP-TSD compared well with experimental data for subsonic 
as well as supersonic freestream Mach numbers which gives 
confidence in the code for aeroelastic prediction.3 


Biiultt anti Pticuiilon 

Results are presented for several configurations to 
demonstrate and evaluate the modifications to the AF 
algorithm of the CAP-TSD code. Calculations are first 
presented for e flat plate airfoil to assess the effectiveness of 
the nonreflecting far field boundary conditions. Calculations 
are next presented for the F-5 wingio and the ONERA M6 
wing ii to demonstrate the improvements due to the 
Engquist-Osher switch, the second-order accurate 
supersonic differencing, and the steady-state convergence 
acceleration. Finally, steady and unsteady results are 
presented for the General Dynamics one-ninth scale F-16C 
aircraft modeli2,i3 to investigate application of the modified 
algorithm to a realistic aircraft configuration. 



( a ) reflecting boundary conditions. 



( b ) nonreflecting boundary conditions. 

Fig. 1 Comparisons of unsteady lift-curve slope for a flat 
plate airfoil at M - 0.85. 
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Elit .Plate. Airfoil -BfliulU 

Unsteady results were obtained for a flat plate airfoil at 
M - 0.85 to test the nonreflecting far field boundary 
conditions. The flat plate airfoil was selected to allow direct 
comparison of results with the exact kernel function method 
of Bland. f 4 The boundary conditions were tested by 
computing the lift coefficient due to the airfoil pitching about 
the quarter chord. Such unsteady forces are typically 
determined by calculating several cycles of forced harmonic 
oscillation with the last cycle providing the estimate of the 
forces. Alternatively, the forces may be obtained indirectly 
from the response due to a smoothly varying exponentially 
shaped pulse. 15 In this procedure, the airfoil is given a 
smalt prescribed pulse in a given mode of motion (in this 
case pitching) and the aerodynamic transients calculated. 


The harmonic response is obtained by a transfer-function 
analysis using fast Fourier transforms. Use of the pulse 
transfer-function technique gives considerable detail in the 
frequency domain with a significant reduction in cost over 
the alternative method of calculating multiple oscillatory 
responses. For the flat plate airfoil, pulse transient 
calculations were performed using 1024 time steps with 
At » 0.2454. The amplitude of the pulse was 0.5o. The grid 
extended 25 chordlengths above and below the airfoil, and 20 
chordlengths upstream and downstream of the airfoil. 
Parallel results were obtained using reflecting (steady- 
state) and nonreflecting far field boundary conditions as 
shown in Fig. 1. The results are plotted as real and 
imaginary components of the unsteady lift-curve slope ck* as 

a function of reduced frequency k. Computations using the 
reflecting boundary conditions, shown in Fig. 1(a), produce 




( a ) steady-state convergence. 


( a ) steady-state convergence. 


NSUP 

NSUP< final) 



NSUP 

NSUP< final) 



( b ) number of supersonic points. 

Fig. 2 Effect of step size on the solution computed using 
the Murman switch for the F-5 wing at M * 0.9 and 
ao - 0o. 


( b ) number of supersonic points. 

Fig. 3 Effect of deleting lime derivatives in the residual on 
the solution computed using the Engquist Osher 
switch for the F-5 wing at M * 0.9 and ao - 0°. 
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oscillations in both the real and imaginary parts for 
0 < k < 0.2. The oscillations are produced by reflected 
disturbances which propagate back into the near field and 
contaminate the solution. When the calculation was repeated 
using the nonreflecting boundary conditions, shown in Fig. 
1(b), the oscillations no longer occur since the boundary 
conditions absorb most of the disturbances that are incident 
on the grid boundaries. Furthermore, these results are in 
excellent agreement with calculations from the kernel 
function method of Ref. 14. 

F-5 Wing Ramils 

Calculations were next performed for the F-5 wing.io to 
assess the algorithm modifications to CAP-TSD. The F-5 
wing has an aspect ratio of 3.16, a leading edge sweep angle 


of 31 .90, and a taper ratio of 0.28. The airfoil section of the 
wing is a modified NACA 65A004.8 airfoil which has a 
drooped nose and is symmetric aft of 40% chord. The F-5 
calculations were performed using a constant step size for a 
total of 500 steps. The freestream Mach number was 
selected as 0.9 and the wing was at 0<> angle of attack. The 
results were obtained to study the steady-state convergence 
characteristics of the modified AF algorithm. The results are 
presented in the form of convergence histories and the 
number of supersonic (NSUP) points versus the iteration 
number. 

In the original AF algorithm of Ref. 4, the Murman type- 
dependent switch was used. Results obtained using the 
unmodified code are presented in Fig. 2. The steady-state 
convergence is shown in Fig. 2(a); the number of supersonic 



( a ) steady-state convergence. 



( a ) steady-state convergence. 




NSUP 

NSUP< final ) 




( b ) number of supersonic points. 

Fig. 4 Effect of over relaxing the residual (with time 
derivatives deleted) on the solution computed using 
the Engquist-Osher switch for the F-5 wing at M - 
0.9 and ao * 0°. 


( b ) number of supersonic points. 

Fig. 5 Effect of deleting all TSD time derivatives on the 
solution computed using the Engquist-Osher switch 
for the F-5 wing at M * 0.9 and ao * 0°. 
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points (NSUP) normalized by the final value are shown in 
Fig. 2(b). For aeroelastic analysis where airloads are 
required rather than pressures, the solution is considered to 
be converged to engineering accuracy when a three to four 
order-of-magnitude reduction in the solution error is 
obtained. The "error*' in the convergence history, as defined 
herein, is the ratio of the maximum |A4| after n iterations to 
the maximum |A4| in the initial solution (first iteration). 
Two sets of results are plotted corresponding to two values of 
step size. At « 0.1 and 0.5. For At - 0.1, the rate of 
convergence is slow and the number of supersonic points 
oscillates about the final value. Increasing the step size to 
At - 0.5 improves the rate of convergence and the 
oscillations in NSUP are significantly damped. The results 
for At * 0.5 also indicate that the number of supersonic 
points is initially more than four and one-half times the 
final value and that "spikes" begin to appear in the 
convergence history after 150 steps. These spikes, which 
represent a numerical instability, are due to a large 
transient caused by the impulsive start from a uniform 
stream using a large step size. If the calculations were 
started with a smaller step size, and then the step size 
increased to the larger value, the numerical instability can 
be avoided. Also, as shown in Refs. 2 and 4, the step size may 
be cycled through very large values such as At - 5.0 to 
achieve faster convergence to steady-state. 

The F-5 calculations with At - 0.5 were then repeated 
with the E-0 switch replacing the Murman switch. These 
results are labeled "unsteady residual" in Fig. 3. The curves 
are identical (within plotting accuracy) to the At ~ 0.5 
curves of Fig. 2 except that the spikes in the convergence 
history are absent. The E-0 switch is more robust than the 
Murman switch and thus the calculation remains stable. 
Furthermore, the rate of convergence to steady-state could 
be increased by deleting the time derivatives in the residual. 
These results, which are labeled "steady residual" in Fig. 3, 
show that after the first 70 steps the solution converges 
faster and the initial overprediction of NSUP is less than that 
computed using the unsteady residual. 

The convergence to steady state could be further 
accelerated by over-relaxing the residual as shown in Fig. 4. 
The results labeled "original residual" are the same as the 
"steady residual" curves presented in Fig. 3. The over- 
relaxed residual results of Fig. 4 were obtained by doubling 
the residual using a - 2.0. These results indicate a faster 
rate of convergence, especially in the first part of the 
calculation, and that NSUP is within 2% of its final value 
after only approximately 50 steps. 

To further investigate the convergence characteristics of 
CAP-TSD, the algorithm was modified to solve the steady TSD 
equation by deleting all of the time derivatives. Calculations 
for the F-5 wing were performed using At - 0,5 and o » 1.0 
to directly compare with parallel results obtained by solving 
the unsteady TSD equation. These comparisons are presented 
in Fig. 5. The convergence history computed using the steady 
algorithm is monotonically decreasing and very smooth in 
comparison with the unsteady algorithm convergence 
history. The steady algorithm solution converges faster and 
does not produce the large initial overprediction of NSUP 
that is characteristic of the unsteady algorithm. The number 
of supersonic points converges rapidly to within 2% of its 
tinal value in only approximately 25 steps. Over-relaxing 
the residual of the steady algorithm also further accelerates 
the convergence to steady-state (not shown). 


QNERA MS Wing Rwulla 

To test the accuracy of the modified CAP-TSD algorithm, 
calculations were performed for the ONERA M6 wing.n The 
M6 wing has an aspect ratio of 3.8, a leading edge sweep 

angle of 30°. and a taper ratio of 0.562. The airfoil section 
of the wing is the ONERA "D" airfoil which is a 10% 
maximum thickness-to-chord ratio conventional section. 
The freestream Mach number was selected as M « 0.84 and 
the wing was at 3.06° angle of attack. These conditions were 
chosen for comparison with the tabulated experimental 
pressure data of Ref. 11. This rather well-known case is a 
very challenging one, especialty for a TSD code, because of 
the complex double shock wave which occurs on the upper 
surface of the wing. 

Steady-state calculations were performed for the M6 
wing by using the AF algorithm with the E-0 switch. The 
results were obtained by cycling the step size through values 
as large as At - 2.0 for a total of 500 steps. This relatively 
large step size corresponds to two root chords of travel per 
time step. A comparison of the resulting CAP-TSD pressures 
with the experimental pressure data is given in Fig. 6 for 
two chords along the span. Results for q - 0.44 are shown in 
Fig. 6(a); results for r\ * 0.65 are shown in Fig. 6(b). The 
data indicate that there is a relatively weak highly-swept 
supersonic-to supersonic shock wave which forms forward 
near the leading edge. The primary supersonic-to-subsonic 
shock which occurs in the midchord region of the wing, 
coalesces with the first shock. Outboard toward the tip, the 
two shocks merge to form a single supersonic-to-subsonic 
shock wave. The CAP-TSD results, obtained using first- 
order-accurate differencing in supersonic regions, are in 
fairly good agreement with the data in predicting the overall 
pressure levels, although differences occur in the regions of 
the shocks. In general, the leading edge suction peak is well 
predicted but the supersonic-to-supersonic shock is 
smeared. When the calculation was repeated using the 
second-order-accurate spatial differencing, a significant 
improvement was obtained in the accuracy of the results. 
The comparisons in Fig. 6 show that the supersonic-to- 
supersonic shock is much more sharply captured by the 
second-order method and consequently the calculated 
pressures are now in very good agreement with the 
experimental data. Calculations were also performed for the 
M6 wing using the original algorithm with the Murman 
switch. These calculations were unsuccessful because of a 
numerical instability which was produced by the highly 
expanded flow about the leading edge of the wing. 

An unsteady calculation was also performed for the M6 
wing at M - 0 84, to investigate the robustness of the 
modified algorithm tor time dependent applications. In this 
demonstration calculation, the wing was forced to oscillate in 
pitch about a line perpendicular to the root at the root 
midchord. The amplitude of the motion was 2° peak -to peak 
about the mean angle of attack of nto * 3.06°. The reduced 
frequency was selected as k * 0.1 and only 300 steps per 
cycle of motion were used. This corresponds to a step size of 
At * 0.1047. Three cycles of motion were computed to 
obtain a periodic solution. Unsteady pressure distributions, 
obtained using first-order and second order accurate 
supersonic differencing, are shown at the maximum pitch 
angle (a - 4.06°) in Fig. 7. Results for rj « 0.44 are shown 
in Fig. 7(a); results for r\ • 0.65 are shown in Fig. 7(b). 
Simitar to the steady-state results, these pressure 
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comparisons illustrate that the supersonic-to-supersonic 
shock is more sharply captured by the second-order method. 
Further instantaneous pressure distributions at two points 
during the third cycle of motion are shown in Fig. 6 for five 
span stations along the wing. Pressures at the wing 
maximum angle of attack (a - 4.06®) and pressures at the 
wing minimum angle of attack (a - 2.06°) are both 
presented in the figure. As the wing pitches up, the shocks 
move aft and the supersonic-to-subsonic shock grows in 
strength. As the wing pitches down, the shocks move forward 
and the supersonic-to-supersonic shock is more sharply 
defined. For this case, both of the shocks oscillate over 
approximately 10% of the chord during a cycle of motion. 
Also, the supersonic-to-supersonic shock at n - 0.80 
periodically appears and disappears during a cycle of motion. 
The results illustrate the large shock motions that the 



x/c 


modified AF algorithm is capable of computing. The 
improved algorithm captures the shocks sharply and is 
sufficiently robust to compute this complex unsteady flow 
using only 300 steps per cycle of motion. 

Ganaral Dvnmlei F-16C Aircraft Modal Remit! 

Results were also obtained for the General Dynamics 
F-16C aircraft model’ 2 to investigate application of the 
modified algorithm to a realistic aircraft configuration. 
Shown in Fig. 9 are the F-16C components that are modeled 
using CAP-TSD. The F-16C is modeled using four lifting 
surfaces and two bodies. The lifting surfaces include: (1) 
the wing with leading and trailing edge control surfaces, (2) 
the launcher, (3) a highly-swept strake, aft strake, and 
shelf surface, and (4) the horizontal tail. The bodies 
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(a) n * 0.44. 



(a) ti - 0.44. 



(b) n -0.65. 

Fig. 6 Effects of first-order and second-order accurate 
supersonic differencing on the steady pressure 
distributions of the ONERA M6 wing at M - 0.84 
and ao - 3.06°. 


(b) r\ - 0.65. 

Fig 7 Effects of first-order and second order accurate 
supersonic differencing on the unsteady pressure 
distributions of the ONERA M6 wing during the 
third cycle of rigid pitching at M « 0.84, ao * 
3.060, a1 „ i.o°, and k - 0.1. 
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Fig. 8 Instantaneous pressure distributions on the ONERA M6 wing during the third cycle of 
rigid pitching at M - 0.84, ao - 3.06°, at » 1.0°, and k « O.t. 


include: (1) the tip missile, and (2) the fuselage. Other 
salient features of the F-16C modeling include 3° linear 
twist washout for the wing, a leading edge control surface 
hinge line that is straight but not of constant-percent chord, 
and 10° anhedral for the horizontal tail. In these 
calculations, the freestream Mach number was M - 0.9 and 
the F-16C aircraft was at 2.38<> angle of attack. Also, the 
leading edge control surface of the wing was deflected 
upwards 2<> for comparison with the experimental steady 
pressure data of Ref. 13. Furthermore, the calculations 
were performed on a grid which conforms to the leading and 
trailing edges of the lifting surfaces and contains 324,000 
points. Since the grid is Cartesian, it was relatively easy to 
generate, even for such a complex configuration as the 
F-16C aircraft. Also, the calculations required only about 
0.88 CPU seconds per time step and thirteen million words 
of memory on the CDC VPS-32 computer at NASA Langley 
Research Center. 

Steady-state calculations were performed for the F-16C 
aircraft using the AF algorithm with the E-0 and Murman 
switches. The E-0 results were obtained using both the 
first-order and second-order accurate supersonic 
differencing. Steady pressure comparisons are given in Fig. 
10 for three span stations of the wing and one span station of 



Fig. 9 CAP-TSD modeling of the General Dynamics one- 
ninth scale F-16C aircraft model. 


the tail. Both sets of E-O results are presented for 
comparison with the experimental data. The results obtained 
using the Murman switch were originally published in Ref. 
2. These results are identical to plotting accuracy with the 
first-order E-O results, and therefore are not shown. The 
steady pressure comparisons indicate that there is a 
moderately strong shock wave on the upper surface of the 
wing and the CAP-TSD pressures agree well with the 
experimental pressures. For the tail, the flow is 
predominantly subcritical and the calculated results again 
agree well with the data. Comparison of pressures computed 
using first-order and second-order accurate supersonic 
differencing shows very small differences. The largest 
difference, for example, occurs on the wing at ij* - 079 
where the second-order calculation predicts a slightly 
stronger shock. 

Unsteady results were also obtained for the F-16C 
aircraft to investigate the robustness of the modified 
algorithm for realistic-aircraft time-dependent 
applications. For simplicity, the calculation was performed 
for a rigid pitching motion where the entire aircraft was 
forced to oscillate about the model moment reference axis at 
a reduced frequency of k * 0.1. The oscillation amplitude 
was chosen as at - 1.5° which is three times the value used 
to obtain similar results presented in Ref. 2. Three cycles of 
motion were computed using 300 steps per cycle of motion 
corresponding to At - 0.1047. Calculations were performed 
using both the Murman and E-0 switches. The solution using 
the original algorithm with the Murman switch, however, 
was numerically unstable for this case as shown in Fig. 1 1 . 
Instantaneous pressure distributions at time steps 94 and 
95 are plotted in the figure, computed using the Murman 
(Fig. 11(a)) and E-0 (Fig. 11(b)) switches. The numerical 
instability begins in the region of the launcher/tip-missile 
where the grid spacing is smallest. Figure 11(a) shows the 
instability in the form of an oscillation in the wing upper 
surface pressure distribution at ij* - 0.94 from 
approximately 30% to 60% chord. The program 
subsequently failed during step 96 which is 21 steps after 
the maximum pitch angle in the first cycle of motion. The 
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Fig. 10 Comparisons between CAP-TSD steady pressure distributions computed using first-order 
and second-order accurate supersonic differencing with the experimental pressure data 
for the wing and tail of the F-16C aircraft model at M « 0.9 and ao - 2.38°. 


calculation involving the modified algorithm (E-0 switch 
with the first-order accurate supersonic differencing) is 
stable, however, as shown in Fig. 11(b). Here the pressure 
distributions for steps 94 and 95 are very similar and the 
calculation proceeds with no difficulty. In fact, the modified 
AF algorithm with the E-0 switch is numerically stable for 
this case with either the first-order or second-order 
supersonic differencing. 

Unsteady pressure distributions along the wing and tail 
during the third cycle of motion are shown in Fig. 12, 
computed using the E-0 switch with the second-order 
accurate supersonic differencing. Two sets of calculated 
pressures are presented corresponding to the aircraft at the 


maximum (a - 3.88°) and minimum (a - 0.88°) pitch 
angles. Comparison of the results indicates that large 
changes in pressure occur along the upper and lower 
surfaces of the wing as the aircraft oscillates in pitch. For 
example, the shock on the wing upper surface oscillates over 
more than 10% of the chord during a cycle of motion. Also, 
the shock is approximately twice as strong at the maximum 
pitch angle as it is at the minimum pitch angle. For the tail, 
the changes in the pressure distributions due to aircraft 
pitching are relatively very small in comparison with the 
changes in wing pressures, as further shown in Fig. 12. The 
tail is located considerably aft of the pitch axis and thus its 
motion is plunge dominated which results in much smaller 
airloads for the low value of k considered. 
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(a) numerical 
switch. 


instability with Murman 


( b ) improved numerical 
Engquist-Osher switch. 


stability with 


Fig. 11 Effect of type-dependent switch on numerical stability for rigid pitching of the F-16C 
aircraft model at M * 0.9, ao - 2.38°. at - 1.5°, and k - 0.1. 
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Fig. 12 Instantaneous pressure distributions on the wing and tail of the F-16C aircraft model 
during the third cycle of rigid pitching at M - 0.9, ct 0 = 2.38°, ai * 1.5°, and k * 0.1. 


Concluding Remarks 

Improvements to a time*accurate approximate 
factorization (A F) algorithm have been implemented for 
steady and unsteady transonic analysis of realistic aircraft 
configurations. These algorithm improvements have been 
made to the CAP-TSD (Computational &eroelasticity 
Erogram * Jransonic Small disturbance) code developed 
recently at NASA Langley Research Center. The AF algorithm 
of the CAP-TSD code solves the unsteady transonic small* 
disturbance equation. The paper described recent changes to 
the code which have significantly improved the stability of 
the AF algorithm and the accuracy of the results . The 
algorithm modifications include: an Engquist-Osher (E-O) 
type dependent switch to treat regions of supersonic flow, 
extension of the E-0 switch for second order spatial 
accuracy, nonreflecting far field boundary conditions tor 
unsteady applications, and several modifications to 
accelerate convergence to steady-state. 

Calculations were presented for the F 5 wing and ONERA 
M6 wing which demonstrated applications of the algorithm 
improvements. The results revealed the superior stability 
characteristics and computational efficiency of the E O 
switch. Much larger time steps were possible using the E-0 
switch, even for comparatively difficult cases. For the 
particularly challenging case of the M6 wing at M » 0.84 and 
ao * 3.06°, the AF algorithm with the E-0 switch was found 
to be stable for time steps as large as At « 2.0. This 
relatively large step size corresponds to two root chords of 
travel per time step. Comparisons of resutts obtained using 
first-order and second-order supersonic differencing 
clearly demonstrated the improved accuracy of the second - 
order method. Changes to the AF algorithm for convergence 
acceleration, namely deleting time-derivatives from the 
original unsteady algorithm and over- relaxing the residual, 
resulted in faster rates of convergence to steady state. 
Converged solutions were obtained in only several hundred 
time steps for the F-5 and M6 wings. An unsteady 
calculation for the M6 wing undergoing a rigid pitching 
oscillation demonstrated the robustness of the modified AF 
algorithm. In this calculation, the shocks oscillated over 


approximately 1 0% of the chord and the flow was computed 
using only 300 steps per cycle of motion. This rather 
difficult case could not be computed using the original 
algorithm. 

Calculations were also presented for the General 
Dynamics one-ninth scale F-16C aircraft model to 
demonstrate application of the modified CAP-TSD code to a 
realistic aircraft configuration. The F-16C components that 
were modeled included: the wing with leading and trailing 
edge control surfaces: a highly-swept strake, aft strake, and 
shelf surface; the tip launcher and missile; the horizontal 
tail; and the fuselage. Steady pressure results at M * 0.9 
and ao * 2.38° compared well with the experimental data. 
Unsteady results were presented for the entire F-16C 
aircraft undergoing a rigid pitching motion with a three 
degree peak-lo-peak oscillation amplitude. The calculation 
was a challenging one for the modified algorithm since the 
flow was computed using only 300 steps per cycle of motion. 
In this calculation, the shock on the upper surface of the 
F-16C wing oscillated over more than 10% of the chord 
which further demonstrates the robustness of the modified 
algorithm. Also, similar to the M6 wing example, this case 
could not have been computed using the original algorithm. 
Therefore, the modifications have significantly improved the 
numerical stability of the AF algorithm and the general 
reliability of the CAP-TSD code for realistic aircraft 
applications. 
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